Bose-Einstein Condensation at a Helium Surface 



E. W. Draeger^ and D. M. Ceperley 
Department of Physics and National Center for Supercomputing Applications, University of 

Illinois Urbana- Champaign, 61801 

o 

^ ! Abstract 

Path Integral Monte Carlo was used to calculate the Bose-Einstein condensate 

in 

CN ' fraction at the surface of a helium film at T = 0.77 K, as a function of density. 

Moving from the center of the slab to the surface, the condensate fraction was 
J — ' found to initially increase with decreasing density to a maximum value of 0.9 

before decreasing. Long wavelength density correlations were observed in the 
' static structure factor at the surface of the slab. Finally, a surface dispersion 

' relation was calculated from imaginary-time density-density correlations. 

a 

^ . 

g . I. INTRODUCTION 



It has been suggested [HI that the condensate fraction in the low density region near 
a ^He surface is significantly larger than the value in bulk helium of 0.1 0. Variational 
Monte Carlo (VMC) simulations by Lewart and Pandharipande of small (A^ = 70) "^He 
droplets using a Jastrow one-body (JOB) trial wave function give evidence for a condensate 
fraction which approaches unity as the density goes to zero in the helium surface. However, 
subsequent calculations by Galli and Reatto have shown that the condensate fraction 
throughout a helium surface computed using VMC is highly sensitive to the choice of trial 
wave function. They found that calculations performed using a shadow wave function with 
a glue term (G-SWF) have enhanced density-density correlations at long wavelengths [§, 
and a maximum condensate fraction of only 0.5. Significant ripplon excitations are expected 
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to reduce the condensate fraction at the surface, as fewer atoms are able to occupy the zero- 
momentum state. Quantum evaporation experiments ||^ can be interpreted as providing 
evidence of an enhanced condensate fraction. 

To avoid the problem of trial function bias and to include finite-temperature effects, we 
have used Path Integral Monte Carlo (PIMC) to calculate the density-density correlation 
functions and condensate fraction at the surface of liquid ^He. We have also used imaginary- 
time correlation functions to calculate the dispersion relation of surface excitations in a free 
helium surface, and find good agreement with experimental thin film measurements. 

II. SIMULATION DETAILS 

Our simulation system consisted of ^He atoms, interacting pair-wise with a very accurate 
potential . Periodic boundary conditions were used with a box size and initial conditions 
chosen to favor a double-sided film oriented perpendicular to the z-axis. To maintain a 
stable film and minimize finite size effects, we added an external potential determined from 
the long-ranged part of the interaction potential and the missing atoms from the other side 
of the slab, so that atoms on each of the two surfaces saw a potential as if they were at the 
surface of a semi-infinite slab . 

Most of the calculations were performed with T = 0.77 K, with an imaginary time step of 
r = 1/20 K^^. We performed simulations of helium slabs containing = 216 and = 432 
atoms, with dimensions 24 x 24 x 17 A and 34 x 34 x 17 A, respectively. Lower temperature 
calculations were also done to determine the temperature dependence. 

III. THE STATIC STRUCTURE FACTOR 

In order to determine the extent to which ripplons are present in a free helium surface, 
we estimated density-density correlations at the surface with the static structure factor 
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defined as: 

S{k\\;z,z') = {pk^^{z)pk^^{z')) (1) 

where pk(z) = J- — Yl G^^'^'^i^i — z), N{z) is the number of particles in the bin at position 
z and k\\ is the wave vector parallel to the surface. This measures the correlation between 
density fluctuations at vertical positions z and z' . 

Fig. shows S{k\\, z, z) curves as a function of density p{z). Each curve is the average of 
both sides of the slab of eight identical simulations (a total of 8 x 10^ Monte Carlo passes). 
At T = 0.77 K, for densities of p{z) = 0.015 and below there is a small enhancement 
of long-wavelength density-density correlations at k\\ = 2'k/L = 0.18 A^^, which is evidence 
for ripplons. However, the curves are closer to the VMC calculations of Galli and Reatto 
which used a JOB trial wave function than those which used the G-SWF form. This does 
not imply that the JOB trial wave function is well-suited to representing an inhomogeneous 
helium system such as a helium slab, but rather the degree to which the G-SWF significantly 
overestimates the effect of ripplons in a free helium surface. Calculations at iV = 216 and 
T = 0.77 K agree with the = 432, T = 0.77 K results of Fig. |^ within statistical error, 
indicating that finite-size effects are negligible. At T = 0.38 K and A^ = 216, we find a 
measurable decrease in the long-wavelength correlations at /cn = 2-k/L = 0.26 A~^. Further 
studies are needed to establish the temperature dependence. 



IV. THE EXCITATION SPECTRUM 

The excitation spectrum can be estimated with path integrals using imaginary-time 

correlation functions 0. The dynamic structure factor is related to the imaginary-time 
density-density correlation function by: 

oo 

F(k,t) = j duoe-'^ S%uj) (2) 

— oo 

= ^ (Pk(t)Pk(O)) . (3) 



To select out the excitations at the free surface, we want to calculate the imaginary-time 
correlation function of propagating surface modes. W. F. Saam |Ty] proposed that the lowest 
quantized hydrodynamic mode (capillary wave) at a free helium surface will have the form 

0fco(r,t) = 0fco(^)e*'^-'-"e-*-'=o* (4) 

where 

0fco(^) oc e-'^W^ (5) 
and the decay constant K,{k) is defined as 

K{k) = -b, + {e + bif' (6) 



h - ^ (7) 



where cxo is the zero-temperature surface tension, po is the bulk density, and s is the zero- 
temperature sound velocity. To calculate the dispersion relations for excitations of this form, 
we use pk in Eq. (H), defined as 



Pk 



J2 e'*^""'" M^i) (8) 



where 



M^) = { (9) 
1 if z > Zc 

and K{k) is defined by Eq. (^). We defined Zc as the point in the surface at which the average 
density was equal to 10% of the bulk density. Tests have shown that the results are not 
sensitive to the value of Zc- 

Extracting the dynamic structure factor by inverting Eq. (Q) is ill-conditioned in the 
presence of statistical noise. It has been shown [|TT|] that a maximum entropy method greatly 
increases the numerical stability. In the maximum entropy method, S{k,uj) is calculated by 
minimizing the function 

-(1/2)Q(S) aCiS) 

HS,a) = '-- x'-—- (10) 



where Q{S) is the "hkehhood" of the PIMC data given an 5* and C{S) is the entropy of a 
given S{k, u) defined with respect to some default model, with a an adjustable parameter 
(also optimized). 

The dispersion energy for a given value of k can be determined from the position of the 
maximum value of S{k,uj). Boninsegni and Ceperley JT2| found that the position of the 



main S{k,uj) peaks for liquid helium agree quite well with experiment, despite significant 
broadening of the excitation spectrum caused by the maximum entropy procedure. The 
dispersion energy of the surface excitations as estimated using this procedure are shown in 
Fig. ^. The two lowest data points, at A; < 0.5 A, had significant fitting error and are only 
qualitatively reliable. Otherwise, we see excellent agreement with the experimental thin film 



data of Lauter et al. [13], including the curvature of the ripplon branch toward the roton 



minimum, proposed as evidence for roton-ripplon hybridization |T5|,|T6 



V. THE BOSE-EINSTEIN CONDENSATE FRACTION 

We define the condensate fraction in the slab geometry by the fraction of atoms at a given 
value of z having precisely A;|| = 0. (Because [k\\,z\ = 0, we can measure the momentum 
parallel to the surface simultaneously with the z-position.) The momentum distribution at 
a distance Zq from the center of the slab is given by 

^^\\i^o) = j dr\\ e^^'^ii-'ii n(r|i;zo) (11) 

where the off-diagonal single particle density matrix is: 

n{Y]z) = \ I dYi---dYN 

p{z)Z J 

X p(ri, r2, ■ ■ ■ Tat, ri + r, r2, ■ ■ ■ tn] P), (12) 

where p is the many-body density matrix and Z = Tr(p) is the partition function. This 
function can be calculated from PIMC []T^ by performing simulations with a single open 



path. We fix the endpoints of the open path at z = zq, and calculate the distribution of 



end-to-end distance n{r\\; zq)- The condensate fraction at a given point in the surface is: 

n{r\\ ^ oc; zo) , . 

M^o) = —I ^ ^- (13) 

Using PIMC, we calculated n{r\\, z) throughout the slab. Nonlinear least-squares fitting 
was used in the region of r\\ < 1.5 A to get an estimate of n{r\\ = 0, z), and n{r\\ cxd, z) 
was calculated by averaging over the region at large r\\ where n{r\\,z) is fiat. At the lowest 
densities, it is not clear whether n{r) has reached its asymptotic limit within the finite 
simulation box. Thus, for = 216 at densities below p{z) = 0.001 A~'^ our results are 
upper bounds to the condensate fraction. 

The condensate fraction no{z) is plotted as a function of average density p{z) in Fig. ^, for 
both A^ = 216 and A^ = 432 helium slabs. As one moves through the surface, the condensate 
fraction initially increases with decreasing density, due to the decreased zero-point motion 
from helium-helium interactions, reaching a maximum value of 0.93(3) at p{z) = 0.002 A~^. 
As the average density decreases below this point, the condensate fraction begins to decrease, 
further evidence for correlated density fluctuations due to ripplons at the surface. This is 
in qualitative agreement with the G-SWF VMC calculations of Galli and Reatto. However, 
the G-SWF trial wave function significantly overestimates the degree to which ripplons are 
present in the surface. 

We find that the probability of a given atom belonging to a long exchange cycle does not 
change if that atom is located above the surface. Hence the density fluctuations pulling atoms 
above the surface are accompanied by other (exchanging) atoms, in qualitative agreement 
with the ripplon model. But there is also present an appreciable fraction (^ 30%) of non- 
exchanging atoms, possibly a result of finite temperature excitations. 

We have presented PIMC calculations of the density- density correlation functions and 
condensate fraction at a free helium surface. These results support the model of a free 
helium surface with ripplons, in which the condensate fraction reaches a maximum at an 
intermediate density in the liquid- vacuum interface, before decreasing at lower densities. Ex- 
perimental probes of the surface will indeed see an enhanced condensate fraction as proposed 
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by Griffin and Stringari 
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FIGURES 

FIG. 1. PIMC density distribution, for both N = 216 and N = 432, T = 0.77 K semi-infinite 
^He slabs (open and solid circles). The two curves are indistinguishable. The dashed line shows 
the effective helium density felt by atoms at z > 0. 

FIG. 2. S{k^^,z,z) vs. throughout the surface region in an = 432, T = 0.77 K slab, 
calculated with PIMC. 

FIG. 3. Dispersion relation of free surface excitations. The results calculated from PIMC and 
maximum-entropy inversion (filled diamonds) for an = 216, T = 0.77 K '^He slab are compared 
with the experimental thin film data (open circles) of Lauter et al. and the DFT results (solid 
line) of Lastri et al. Also shown is the dispersion relation of bulk ^He at SVP (dashed line). 

FIG. 4. n(r||,2:) vs r||, throughout the surface region of = 216 (top) and A'" = 432 
(bottom) ^He slabs at T = 0.77 K, calculated with PIMC. Density labels correspond to 
z = 0,8, 9, 10, 10.5, 11, 12, 13 A for A^ = 216, and z = 0,9, 11, 12, 13 A for A^ = 432. 

FIG. 5. Condensate fraction vs. density throughout the surface region of a T = 0.77 K '^He 
slab, calculated from the n(r||,z) curves in Fig. ^ Also shown are the VMC calculations of Galli 
and Reatto. 
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